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1 Introduction 



Semiclassical approximations of the solutions to the wide variety of wave equations 
in mathematical physics are well known since at least the nineteenth century, mainly 
based on the ideas of short wavelength approximations. Their applications became 
of extreme importance in solving the energy eigenvalue problem of the stationary 
Schrodinger equation. In fact, this has been done by the famous EBK quantiza- 
tion, initiated by Bohr and Sommerfeld for the separable potentials, generalized by 
Einstein to the integrable potentials (by quantizing the classical invariant tori), and 
completed by Maslov (by taking into account the phase corrections at the classical 
turning points and caustics in more than one degree of freedom). This is known 
under the name torus quantization. Einstein's work was in fact completed in 1917 
even before the quantum mechanics was discovered, and it has been later realized 
that it is (after Maslov phase corrections are taken into account) precisely the lead- 
ing order semiclassical approximation to the exact solution. Einstein knew from the 
communications with Henri Poincare that the integrability in classical Hamiltonian 
systems is exceptional, and thus that nonintegrability and nonexistence of invariant 
tori is the generic rule. How to proceed in such case was a mystery for him. It is 
later through the works of Gutzwiller around 1970 (see his book Gutzwiller 1990) 
that we got an answer to this question, which is a systematic approach to represent 
the leading semiclassical approximation as a sum over contributions associated with 
classical periodic orbits, the so called Gutzwiller theory. All these findings were of 
crucial heuristic, qualitative and to large extent also quantitative value and impor- 
tance especially in the course of studying the quantum chaos, which is the study of 
the quantum signatures of classical chaos. 

On the other hand in 1-dim systems more can be done and indeed there is a quite 
well developed theory (see e.g. Delabaere et al (1997), Balian et al (1979), Voros 
(1983) and references therein). Nevertheless if we ask for explicit results not so 
much is known. However there are the systematic semiclassical procedures in the 
one dimensional potentials, which are of course always integrable, and where the 
leading torus approximation is only precisely the first term in a certain WKB series 
expansion, where each term can be obtained recursively, in principle, and where in 
certain cases the series can be summed and it can be shown to yield precisely the 
exact result. In this regard a very important classic and historic paper is by Bender 
et al (1977). See also (Robnik and Salasnich 1997a,b) and Romanovski and Robnik 
(1999) and (Salasnich and Sattin 1997). The Schrodinger eigenvalue problem is very 
difficult even for general one dimensional potentials and any analytically tractable 
WKB approximations and exact results are of extreme importance in qualitative 
and quantitative understanding of general solutions. Also, they provide a basis for 
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the many dimensional integrable cases, where almost nothing is known beyond the 
leading Einstein-Maslov term. In the following we shall use the terminology "WKB 
series expansion" as the synonym of the " semiclassical expansions" which we clearly 
define below. 

Although at present the WKB theory is very deeply developed and its methods 
are very important for many applications, there are only a few works where the 
problem of effective calculation of the terms of WKB expansions is discussed. In this 
direction a pioneering paper is by Bender et al (1977), where the authors investigated 
the structure of the terms of WKB expansions and also applied the methods to 
compute the eigenvalues of the potential V(x) = x N (N even). In our present 
paper we perform further study of the problem of effective computation of WKB 
series started in Bender et al (1977). We obtain new recurrence formulae for WKB 
terms for arbitrary analytic potentials and for the polynomial potential V(x) = 
x N . The known algorithms are very laborious, because they involve operations of 
differentiation and collection of similar terms in polynomials, which are extremely 
time-consuming as the order increases, even when modern computer algebra systems 
are used. In contradiction, in computing by means of our recurrence formulae one 
does only arithmetic operations with rational numbers, as we have explicit formulae 
for the WKB terms, except for the numerical coefficients. 

Using the obtained functional form for the general WKB terms a' k we give a sys- 
tematic algorithm for a dramatic simplification of the integrated WKB terms § a' k dx 
that enter the energy eigenvalue equation, by showing systematically that the idea 
of Bender et al (1977) can be fully implemented resulting in the maximum pos- 
sible simplification of the integrated terms. This is achieved by finding terms of 
the integrand which are complete derivatives of some function and thus give zero 
contribution when integrated round a closed integration contour in the complexified 
coordinate plane x. We also derive almost explicit formulae for the WKB terms for 
the energy eigenvalues of the homogeneous power law potentials V(x) = x N , where 
N is even. 

These results go substantially beyond the results of Bender et al (1977) and indeed it 
should be emphasized that the main algebraic ideas behind our present work are due 
to certain remarkable similarities between our present problems and those involved 
in calculating the normal forms and Lyapunov focus quantities in the power law 
differential equations (of one degree of freedom) and maps (see Romanovski 1993, 
Romanovski and Rauh 1998). There are also some common features between the 
problem of reduction of the coefficients a' k of the WKB series and the problem of 
finding a basis of the ideal of Lyapunov focus quantities (the so-called local 16th 
Hilbert problem, see Romanovski 1996). Some introduction to the WKB method 
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can be found in (Bender and Orszag 1978). 



We consider the two-turning point eigenvalue problem for the one-dimensional Schro- 
dinger equation 

A 2 

dx 



h 2 ^ + V(x)]ij(x) = Eij(x). (1) 



We can always write the wavefunction as 

${x) = exp|^a(a;)| (2) 
where the phase a(x) is a complex function that satisfies the differential equation 

a' 2 (x) + ha"(x) = (V(x) - E) = f Q(x). (3) 
The WKB expansion for the phase is 

CO 

= $>V fc (x). (4) 

k=0 

Substituting into (|J) and comparing like powers of h gives the recursion relation 

i n— 1 

< 2 = Q(x), o f n = -— r (£a' k o f n _ lt + ai. 1 ). (5) 

za o k=l 

Computing few first functions a' k by means of the recurrent formula we get 

/^7T / ~Q'( X 



, 5Q'(xY -AQ(x)Q"(x) 

°2 = —6 • 

32Q(x)* 

15 Q'{xf + 18 Q(x) Q'(x) Q"(x) - 4 Q(xf ( 



(7) 



x 



a ' 3 = R A n( \ A ^ 

§AQ{x) 



and 



a' A = (1105 Q'(x) 4 - 1768Q{x)Q'{x) 2 Q"{x) + AA8Q{x) 2 Q'{x) Q (3 \x) + (9) 

304Q(x) 2 Q"{x) 2 - 64Q(x) 3 Q (4) (x))/2048 Q(x)% . 
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For the analytical potential V(x) the following quantization condition is known (see 
Dunham (1932), Froman and Froman (1977), Fedoryuk (1983)): 




where n q > is an integer number and 7 is a contour surrounding the turning points 
on the real axis. This relation is an equation with respect to E and using it one can 
find the asymptotic of the eigenvalues E n {Ti) (see e.g. Balian et al (1979), Fedoryuk 
(1983) and references therein). In some cases the series (|1(]) can be summed exactly 
(see Bender et al (1977), Robnik and Salasnich (1997a,b), Romanovski and Robnik 
(1999), Salasnich and Sattin (1997)). 

The zero-order term of the WKB expansion is given by 

i j> da = J dx^E-V(x), (11) 

the first odd term is 

= — (12) 

and to find the higher order terms we need to compute the functions a' k using the 
recursion relation fl5|). We note that the odd-order terms (except the first order 
for a[) yield integrals that vanish exactly, because, as it follows from the results of 
Froman (1966), the functions <j' 2k+1 are total derivatives. 



2 An algorithm for computing a' k 

We will look for a general formula for the functions o' k . As it is known one of the most 
powerful tools for investigation of recurrence relations is the method of generating 
functions (see e.g. Graham et al 1994). The most widely used in combinatorics 
generating functions are the ones with a single variable, for example, a generating 
function for the sequence {go, g±, g%, . . .} is the formal series 

G{z) = Y.9nZ n . (13) 

n>0 

For a multi-index sequence {g( yi ,...,y m ))}y i eN, where N is the set of non-negative 
integers, a generating function has the form 

G( Zl , ...,z m )=J2 9(y 1 ,-,y m )4 1 .--z%?. (14) 

Vi>0 
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Moreover, we can also consider a sequence, where every term has only finite num- 
ber of indices, but the total number of indices is unbounded (e.g. 9( yi ), 9{yi,y 2 )i 
g( yi ,y 3 ,..., y „), ■ ■ •) with the generating function 

G(zi,z 2 , • • •) = J2 9^, (15) 

where M = U^ =1 N fc , z = (z\ , . . . , z s ) and z 1 = zj 1 ... z] a . Thus in the last case 
G(z\, Z2, . . .) is an element of the ring of formal power series in the infinite number 
of variables, z\,z 2 , . . . . 

We now apply the method of generating functions to computing of the WKB expan- 
sion for the phase. With any vector v = (1/1,1/2, ... , v\) we associate the operator 

L(y) = \-v 1 + 2-v 2 + ... + l-v h (16) 
where the vector v runs through whole M. 
We will show that the functions a' m are of the form 

(T m= 2^ 3m-l , (17) 

u:L{u)=m 2 

where for a vector v = {u u . . . , v{) we denote QW = (Q') Vl (Q") U2 . . . (Q {l) ) Vl , \v\ = 
V\ + . . . + V\ and the coefficients U v satisfy the recurrence relation 

Ur l UfU> + (4 - L( V ) - n) + g (filig^l , m 

^ fj,,e^o,fi+e=u 4 i=i ^ 

where Uo = — 1 and we put U 1 = if among the coordinates of the vector 7 there 
is negative one, and we denote by u(i) (i = 1, ...,/ — 1) the vector . . . ,Vi + 
1, u i+ i - 1, ... , ui). 

It should be mentioned that the complexity of functions a' n increases rapidly with 
the order n, and it is remarkable that applying our almost explicit formulae (|T7|), 
( |l"8| ) we can go much further than by using the well-known recursion relation (||) 
(see the Appendix). 

Note that the number of solutions v of the equation 

L{v) = m, (19) 

where m > and v is a vector with non-negative coordinates, equals the number 
of partitions p(m) of the integer m (there exists a theory and a formula for p(m); 
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see e.g. Andrews (1976)). Therefore as a corollary of formula (|T7|) we find that the 
number of terms of the function a' k cannot exceed p(m). This fact was for the first 
time observed by Bender et al (1977). 

We prove the formulae (|TTD, ( |TBD by induction on m. Indeed, for m = 1,2 the 
statement holds. Let us suppose that it is true for all m less than k. Then for 

m = k we get 



2a' 



E 

e,fj,:L(d)=i,L(fi)=k-i 



UeUnQt-M-mQW) 
2Q^ ' 



Obviously, 



L(9) = i, L(fj) = k - i L(9 + //) = k, 



(20) 



(21) 



and here we assume that the dimensions of vectors 9 and \x are the same, if not 
we simply extend the dimension of the smaller one by putting zeros for the excess 
coordinates. On the other hand, taking into account that all coefficients of the 
operator ([IB]) are positive, it is easy to verify that 



L(9 + f x) = k ^L(9)=j,L(fx) = k-j, < j < k 
(this property is the crucial one for the presented method). Thus (2~0)- 



k— 1 _./_./ 

E °i°k-i 
9rr' 



E 

0,H:L(0+fi)=k 



2Q 3J ^ 



(22) 
yield 

(23) 



i.e. we get an expression of the form QT7D. 

For the last term of the recurrence formula (||) we get 

"(2 - k - 2| A t|)Q fc -M- 1 Q(Mi+i>M2 > .-,Mfc-i) 



u k-i 
2a' n 



E 



-U„ 



4Q 2 



+ 



Q 



(Q iM) y 



/j,:L(ii)=k-l 

Note that for the vector /i = (fix, . . . , fik-i, 0) 

[Q M ]' = E^Q (A(i)) , 



2g ; 



(24) 
(25) 

where we denote by //(z) {i = 1, . . . , k — 1) the vector (/i!,...,/^ — + 
1, . . . ,/i fc ) (z = 1, . . . , k — 1). It is readily seen that L(fi(i)) = L(/i) + l = k, therefore, 
formulae ©, (H)-(|2|) yield that (0), (0) hold. □ 



i=i 
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Using the recurrence relation (18) we can obtain the differential equation for the 
generating function of the sequence U v 

U(z) = U( Zl ,...)= Y, U »z V - ( 26 ) 
Let us rewrite (118!) in the form 



1 3 

U(u u ...,ui) = r U v u e + ^(i/i-i,^,...,^) - -fV\U{v x -\,v?.,...,v{) ~ (27) 

- 2> + 2)^C/ K _ 1)I/2i ...^ ) + ^ ^ - [i/ = 0], 

4 i=2 i=l Z 

where [a = (3] denotes the function, which equals 1 if a = (3 and otherwise. 

Using obvious properties of generating functions (see e.g. Graham et al 1994) we 
get from (|27|) 



U = \{U + lf + Zl U - \zMU)' zl - £ l -^z xZi U' Zi + \ £ - 1. (28) 

It means if we fix any integer I and, therefore, the variables z±,...,zi, then the 
function 

U(z 1 ,...,z l ) = U(z 1 ,...,z l ,0,0,---) (29) 
is the solution of the equation (p8|) with the initial conditions 

tf(0) = -1,^,(0) = -^ (30) 

(we get the initial conditions from ( |i~8| ) taking into account that C/(o,...,o,i) = ^ or 
the vectors with the only ith coordinate different from zero). So, the coefficients U u 
that enter the functions a' k ( |T7| ) are precisely the coefficients of the Taylor expansion 
of the function U defined by the differential equation (|28|) with the initial conditions 

Coefficients of the form U(n,o 0) depend on the coefficients of the same form. There- 
fore we get from (|28| ) that the function U(z) = J2"^Lo U n z n satisfies the differential 
equation 

U = 1 -{U + lf + zU - \z{zU)' z -l, (31) 
which is the Riccatti equation 

3z 2 U' z - 2U 2 - zU + 2 = 0. (32) 
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Note that as an immediate corollary of formula (|T7|) we get that for the harmonic 
oscillator, i.e. when Q = x 2 — E, the WKB series ( |TUD terminates after the first two 
terms, namely, 

/ da k = (33) 
for all k > 2. Indeed, in this case (|17| ) yields 

[m/2] jj cym~i ^m—H 

, ^ V(m-2i,i)* X 

U m i-^ =:3m-l-2i ' \°^> 

4=0 y/X 2 - E 

where [m/2] stands for the integer part of m/2. 
It is obvious 

Res^ — Zm _ x _ 2l = (35) 

Vr — E 

for all m > 1. Therefore (^) takes place. 

3 An algorithm for the simplification of the func- 
tions dak 

It was pointed out by Bender et al (1977) that the functions dcr^ can be dramatically 
simplified by adding and subtracting total derivatives (obviously, such operation 
does not change the integrals (|T0|) ). They also carried out numerical experiments to 
obtain different simplifications of these functions. It is easily seen that the formula 
( |r?D provides an effective way to reduce the number of terms in the expressions for 
the functions dak- 

We will look for a function of the form 

H-.L(fi)=k Q 2 

where are to be determined, such that 

|n = <+. (37) 

By comparing coefficients of ^ + 1^+?^ ■ , in both parts of (|37D we get the system 

3 ~ L(Z l ) ~ 2|t/| W^-i,^,...^ ) + + 1)W„(0 = E^, (38) 

2 4=1 
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where v — . . . , v{) runs through the whole set of solutions of equation ( |19D and 
Up are defined by the recurrence relation (|TTD , (|T8"D . 

Thus to simplify the function a' k+1 one can solve the system fl38|) of p{k+l) equations 
in p(k) variables W^. For example, to simplify a' A we write down the corresponding 
system (|38D and get 

U 4 = -§W 3 

t/( 2) i) = 3W 3 - lW {1>1) 

U{i,o,i) = W {1<1) - |W (0 ,o,i) (39) 

^(0,2) = 

^(0,0,0,1) = ^(0,0,1) 

where U 4 = 1105/2048, C/ (2)1) = -1768/2048, £/ (0 , 2 ) = 304/2048, £/ (li0)1) = 448/2048, 
^(0,0,0,1) = —64/2048. We see that the matrix, corresponding to the first three equa- 
tions is the triangular one, so we can kill three terms in the expression for da A . 
Computing we get 

1105 991 93 

W 3 = , Wn i) = — , W, 001) = -—. (40) 

3 9216 (,) 1536 { "> 768 V ; 

Hgiicg 

; d _ 7Q"(xf-2Q(x)Q^(x) 
° A dx 3 1536Q(a;) 7 /2 ' 1 j 

in accordance with Bender et al (1977). 

Let us denote by p(k) the number of partitions of k which contain at least one 
1. Obviously, p(k) = p(k + 1). It appears that the optimal strategy to simplify 
a' k is as follows. In system ( |38| ) where L(v) = k we consider the equations with 
v such that v x ^ 0. There are p{k) = p{k — 1) such equations and according to 
(|36D we have exactly pik — 1) variables. It turns out that we can always write the 
systems with U v such that ^ / in the triangular form (like system (^)). To 
see this we set the following order on vectors of N': we say that v < \x if the first 
nonzero entry from the left in \x — v is positive (this order is known in computational 
algebra as the lexicographic one) . Then if we write down the equations of the system 
(|38|), corresponding to U u a), U u <2), ... in the decreasing order z/W > u^ 2 \ . . . and the 
variables in these equations also in the decreasing order then we get that the 
matrix corresponding to the first p(k — l) equations is the triangular p(k — 1) xp(k— 1) 
matrix (because v(i) > (v\ — 1, z/ 2 , . . . , vj) for 1 < i < I — 1). Moreover, the 
diagonal elements of the matrix are equal to (3 — L{y) — 2\v\)/2 and, therefore, are 
different from zero. Hence, we get that after the simplification a' 2n contains at most 
p(2n) — p(2n — 1) terms. 
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To end this section we show that in some cases we can replace calculation of contour 
integral by computing a Riemann integral, namely we shall show that formula (ff3|) 
below applies. 



First we note that taking into account that Q = V(x) — E we can write for even m 
formula fllTD in the form 



v:L(y)=m 



[m 



3 + 2\u\)\\^Ef- 1+ ^^/E~^v' 



Let us now suppose that Q = V(x) — E, where V(x) is an analytic function with 
one minimum and V'(x) 7^ 0, if x 7^ 0. We will show that 

d(T ™ = 2 E , o|„ m , a ^-l+IH / ,/^^7 ^ ( 43 ) 



v:L(y)=m 



im 



3 + 2\v\)\\ dEf-^+M J X1 ^E~^V 



where V(x\) = V(x 2 ) = E , x\ < x 2 . Taking into account fl42|) it is easy to see that 
to prove ( f4"3|) it is sufficient to show that 

,dM^7E^ dx = 2 dM^L VE=v dx ' (44) 

Note that due to the theorem on differentiation upon a parameter (see e.g. Sidorov 
et al 1976) if 

F(E) — <f f(x, E)dx (45) 

and 

1) 7 is a finite piecewise- smooth curve; 

2) the function f(x, E) is continuous with respect to (x, E) for x G 7, E G £>, where 
-D is a domain of the complex plane; 

3) for every fixed 167 the function f(x, E) is analytic upon E in D, then F(E) is 
analytic in D and 

for EeD. 



Let us cut the complex plane between the turning points 07 and x 2 to get a single- 
valued function and fix the contour (x\ + p, x 2 — p) U <7 U (X2 — p, 27 + p) U C2, where 
p is small ci, C2 are circles : c\ = x 2 + pe 1 *, c 2 = x\ + pe 1 * and t G [0, 27r]. Then the 
conditions l)-3) are satisfied (with D being a small neighborhood of E ). Therefore 

«9f-i+M af-i+kl /■ yM 

7 9£?t -i+IH y/E^V** = dEf~^\ ly VE^V dX = (47) 
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Let us denote 



(2 



x-o-p y(W y yM y y"M 

-xo+p VE-V Lx VE-V Jc2 VE-V 



dx). 



9(E) = I 

Jci 



-dx. 



, (48) 

Due to the theorem mentioned above g(E) is analytic in a neighborhood of Eq, 
therefore 

Noting that 

y V {u) f 
\g(Eo) \ = f —j===ydx = f ; dx < const p 



g(E)^g(E )+g\E )(E-E ). 



VEV^v 



c l yJV'(Xi)(x-X 1 )+. 



1/2 



(49) 



(50) 



we obtain that p(-E') — > when p —>■ 0, E —>■ E . It means that formula ([14]) 
indeed holds. This formula has been found useful for computing the WKB series for 
Coulomb potential and the potential V(x) = Uo/cos 2 (ax) (Robnik and Salasnich 
1997a,b). 



4 Potentials of the form V(x) = x 



,N 



In recent decades many studies have been devoted to the investigation of the semi- 
classical expansions for the potentials of the form V(x) = x N and important results 
have been achieved (see e.g. Balian et al (1979), Voros (1983) and references therein). 

One of the basic formulae for these potentials was obtained by Bender et al (1977) 
and is as follows: 

-r oo 

vrK + ~) = E 1/N+1/2 £ £^ (1+2/A °a„(A0, (51) 

^ n=0 

where N is an even integer number, 

a (N) - (-ir^-v^rq + ^p^) 



(2 + 2n)\T(^ + 



N 



and P n (N) are polynomials of the variable A" with integer coefficients. The first eight 
polynomials P n {N) where computed using MACSYMA computer algebra program 
by Bender et al (1977) for the general potential V(x) = x N and, as it was mentioned, 
the computation of the eighth polynomial has already faced difficulties. When A" 
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is a fixed integer number then instead of the polynomials a n we have integers. The 
special case of quartic oscillator (N = 4) has been deeply investigated by Balian, 
Parisi, Voros and others. In this case the expression of the type fl5"T|) is written 
(Balian et al 1979) in the form 

1 oo 

27r(n q + -)h=J2b n cT l ~ 2n h 2n , (53) 

Z n=0 

where 

a = £ 3/4 5(3/2, 1/4) (54) 

is the classical action around the closed orbit of energy E (here B(x,y) is the beta 
function) and b n in this case are rational numbers. As it is reported by Balian et 
al (1979) with REDUCE language they were able to compute the first seventeen 
coefficients b n . Then they had to switch to the ordinary numerical computation and 
computed b n up to n = 53 (in (Voros 1983) the results of computation up to n = 60 
are presented). 

It was mentioned in Balian et al (1979) and in Bender et al (1977) that the authors 
do not know any closed form or a simple law for the coefficients a n (N) and b n . In 
this section we partially answer this question. Although we also were not able to find 
any closed form for the functions a n (or numbers b n ) we obtain a simple recurrence 
formula, where only operations of summation and multiplication of rational numbers 
are involved, whereas using the usual way of the above mentioned papers one needs 
first to compute rational functions of the form f(x, \JE — x N ) and then evaluate 
contour integrals. We carried out computer experiments and found out that using 
our algorithm with Mathematica 4.0 on a PC with 128 MB RAM we were able to 
compute in closed arithmetic form the coefficients b n at least up to n = 190 for the 
quartic potential (see also the Appendix). 

It can be proven (Robnik and Romanovski 2000) by induction using the recursion 
relation (|) that for the potential V(x) = x N the coefficients a' k {k > 1) of the WKB 
expansion have the form 

( -\3k-l -k+N k-1 

°i = E -U , ,,/^ '<"' V . (55) 

j=0 



(E - x N ' 



where we choose sjE — x n = i\/x n — E, and the coefficients Ak-j-ij of the mono- 
mials E k ~i~ 1 xi N are computed according to the recurrence formula 

1 s l-l 

A s,i = ~ E E A i>j A s _ i>l _ 1 _ j + (56) 
l(2 + N) + (2 + 3N)s- N (N-l)l + N-s. 

A^_i H Ag-lJ. 
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with A 0fi = N/4 and 

A a3 = if a < or (3 < 0. (57) 

Using equation (^) we can get the differential equation for the generating function 
of the coefficients A S) i, and in the special case A SjQ one can find the explicit formula 
(Robnik and Romanovski 2000) 

iV> 

As >° = 2^(N-s-l)V (58) 



For even k, (k <-> 2k) we can write formula ( 55]) in the form 

2k-l T)3k pfik 

«4 = E jjppw^^""^^ - -"I" 2 - (59) 

As above we can replace the integration on a contour with the integration between 
turning points, then, taking into account that 

x a -\a e -x e f- 1 dx = — — B(~,/3), (60) 

where a, 9, Re a, Re (3 > 0, and B is the beta-function, and noting that 

r(f + , + ^ + i) 



T(l + s + ^-3k) 

,3 1-2/^,3 l-2k 1N ,3 l-2k OI , 

^ + S + ^v- )( 2 + S + ^v-- 1) --- ( 2 +S + nv-- 3 ' ) (61) 
we get from fl5"5| ) 

da 2k = 2 da 2k = 

l-y J-E 1 /™ 

(6A;-3)!!iV E2 g ^ r(| + s + 1 - 3k + ^) ' (62) 

and using the equality T(l + z) = ^r(z) we obtain finally the coefficients of the 
WKB expansion 

— > + E nc + ^) 2 if '(^ + ^ - (63) 

i=l s=l s=l 

where k > 1 and y4 2 fc~j-i,j are computed according to (|56"D and 

2i£*+* v ^r(i + i 
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5 Conclusions 



In this paper we have investigated the WKB approximations as a series for arbitrary 
analytic potentials. Following the classic paper by Bender et al (1977) we have 
achieved the following results: (i) We find the explicit functional form for the general 
WKB terms a' k , where one has only to solve a general recursion relation for the 
numerical rational coefficients, (ii) We give a systematic algorithm for a dramatic 
simplification of the integrated WKB terms § a' k dx that enter the energy eigenvalue 
equation, (iii) We derive almost explicit formulae for the WKB terms for the energy 
eigenvalues of the homogeneous power law potentials V(x) = x , where N is even. 
In particular, we obtain effective algorithms to compute and reduce the terms of 
these series. In computing by means of these formulae we manipulate only with 
numbers and do not need to collect similar terms of a polynomial, which we must do 
otherwise when we use just the recursion formula (§[). Application of the obtained 
formulae along with the reduction formula fl38|) considerably simplifies calculations, 
especially if we need to compute high order terms. 
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Appendix 

Here we present the results of computer experiments which we carried out with 
Mathematica 4.0 on our PC with 450 MHz processor and 128 MB RAM to compare 
the efficiency of our algorithms based on equations (|T7|), fll8|) and fl56|) , (|63|) with 
traditional ones. 

To compute a' n by Mathematica using formula @ one can just use it in the form 
presented in the text, but for computing by means of formulae fllTD, fll8Tj a procedure 
has to be written. 
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Table 1: (a) the order n of the coefficient a' n , (b) the number of terms in the 
coefficient a' n (equals p(n)), (c) CPU time in seconds for computing a' n by means of 
formula (|5]), (d) CPU time for computing a' n using formulae (|T7D, (|18|). 



(a) 


(b) 


(c) 


(d) 


(a) 


(b) 


(c) 


(d) 


10 


42 


0.47 


0.54 


30 


5604 


18,542 


1,402 


15 


176 


5.7 


5.2 


34 


12310 


114,219 


4,981 


20 


627 


83 


41 


35 


14883 




6,733 


25 


1958 


1,509 


256 


40 


37338 




29,940 



In the case of the potential V(x) = x A for computing § cr' 2n dx one can use formulae 
(|56|), (|63| ) precisely in the form presented in the paper, but it is necessary preliminary 
to define A-\^ = = 0. 

In (Voros 1983) the table of results of numerical calculations of the numbers 

2-1— TT^ CSC(^|^) 

n ~ r(2n-l) 

which should not be confused with our A 8 j and where the numbers b n are defined 
in equation (|53"D, with an estimated accuracy of 34 digits are presented and it is 
mentioned there that the accuracy is not guaranteed for n close to 60. Indeed, we 
found perfect correspondence with the table of (Voros 1983) up to n = 52, however 
for larger n there is a disagreement with our calculations, presented in Table 2 
(digits which differ from those obtained by Voros (1983) are underlined). It should 
be emphasized that here we calculate A n in the exact arithmetic form, which includes 
rational numbers and the gamma function, but show here the numerical results just 
for the purpose of comparing them with those of (Voros 1983). 

We computed the coefficients b n in exact arithmetic form according to formulae Q5TS|), 
fl63|) up to n = 190 and it took 143555 sec CPU time to do so. It was possible to 
continue computations according to the memory capacity, however the computations 
became too time consuming. As it is reported by Balian et al (1979) with REDUCE 
language they were able to go up to n = 16. 
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Table 2: V(x) = x 4 



53 0.9949789006826957939983352522204201 

54 -0.9950733730448721436038332477930891 

55 0.9951643566345461007309335419123495 

56 -0.9952520411721716125492612513022471 

57 0.9953366028692452324220366394 796240 

58 -0.9954182056094612326949059 324394196 

59 0.9954970020080866308729 698916050341 

60 -0.9955731343639551368199525397204894 
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